rm(list = ls())
library(pacman)
p_load(vars,matlab,magrittr,reshape2,stringr,TBHB,foreach,GVARX,ggplot2,tidyverse)
devtools::load_all()

load('./data-raw/regdata2017.rdata') # From ReadyForData2017.R
# from ProEU.r 覆盖了regdata2017.rdata中的unctnMain,增加了欧阳的经济不确定性今天和明天
load('./data-raw/CHNEU.rdata')
unctnMain[,-1] <- apply(unctnMain[,-1], 2, log)
regdata$USA_unctn <- log(regdata$USA_unctn)

load('./data-raw/epuCHN.rdata') # 中国EPU数据
chn$epuCHN <- log(chn$epuCHN)
load('./data-raw/indus_tfp.rdata') # 中国行业TFP
load('./data-raw/ExpImp_weight.RData') # 计算外国变量的外部加权矩阵
load('data-raw/pat.rdata') # 27个国家的专利数据, from ClrPat.R

# from level to TB or HB
TBvar <- names(regdata)[str_detect(names(regdata),'gdp')]
for (i in TBvar) {
  regdata[,i] <- FromLevelToTB(regdata,i,'TB')
}
HBvar <- names(regdata)[str_detect(names(regdata),'cpi')]
for (i in HBvar) {
  regdata[,i] <- FromLevelToTB(regdata,i,'HB')
}

EPUtype <- 'eco'
indus <- ''
patnm <- ''

# ---------注意修改此处------------
# replace EPU: 美国不同政策EPU互换
if (FALSE){
  EPUtype <- tovar <- 'trade'
  regdata <- ReplaceVar(regdata, fromvar = 'USA_unctn', todata = unctnMain, tovar = tovar)
}
# 替换CHN的cpi成EPU for mechanism analysis
if (FALSE){
  regdata <- ReplaceVar(regdata, fromvar = 'CHN_cpi', todata = chn, tovar = 'epuCHN')
}
# 专利数据替换
if (T){
  # 选择部分国家
  sele_cnt <- names(pat)[-1]
  # sele_cnt <- sele_cnt[!(sele_cnt %in% c('SAU'))] # 删国家
  # pat <- pat[, c('date',sele_cnt)]

  sele_cnt_med <- str_detect(names(regdata), paste(sele_cnt,collapse = '|'))
  sele_cnt_med[1] <- TRUE
  regdata <- regdata[,sele_cnt_med]

  patnm <- '_pat'
  for (i in 2:ncol(pat)) {
    regdata <- ReplaceVar(regdata, fromvar = paste(names(pat)[i],'TFP',sep = '_'),
                          todata = pat, tovar = names(pat)[i])
  }
}
# 行业替换
if (FALSE){
  # replace Chinese TFP with industrial TFP
  indus <- tovar <- '造纸及纸制品业'
  regdata <- ReplaceVar(regdata, fromvar = 'CHN_TFP', todata = tfp, tovar = tovar)
}
# ---------------------------------

regdata <- na.omit(regdata)

# 变成国家，时间，gdp,cpi,tfp,unctn格式的data.frame
regdata <- melt(regdata, id.vars = 'date')
regdata$ID <- str_split_fixed(regdata$variable,'_',2)[,1]
regdata$varname <- str_split_fixed(regdata$variable,'_',2)[,2]
regdata <- dcast(regdata, ID + date ~ varname, value.var = 'value')
names(regdata)[2] <- 'Time'
regdata$Time <- str_replace(regdata$Time,'M','/') %>% str_c('/01')

# 外国变量构造
TradeWeight <- TradeWeight[1:17]
for (i in 1:length(TradeWeight)) {
  TradeWeight[[i]] <- TradeWeight[[i]][,-1] %>% as.matrix() %>% t()
  # 选择部分国家
  colnames(TradeWeight[[i]]) <- rownames(TradeWeight[[i]])
  TradeWeight[[i]] <- TradeWeight[[i]][rownames(TradeWeight[[i]]) %in% sele_cnt,colnames(TradeWeight[[i]]) %in% sele_cnt]
}

# TradeWeight <- TradeWeight[2:18]
FV <- GVAR_Ft(regdata,TradeWeight[[5]])
names(FV) <- unique(regdata$ID)

# 设置内生变量，及对应滞后
endovar <- c('cpi','gdp','TFP')

# get W in appendix A.16 of GVAR toolbox. A link matrix.
cntorder <- unique(regdata$ID)
# xt <- NULL
# for (i in 1:length(cntorder)) {
#   xt <- paste(cntorder[i],endovar,sep = '_') %>% c(xt,.)
# }
# W <- GetW(xt,TradeWeight[['2000']])

# 计算GVAR: p is VAR lag, q is foreign variable lag, trim is ratio truncation
p <- 1
q <- 1
countryrlt <-TGVAR(regdata,endovar, thrsholdvar = 'gdp', trim = 0.3, p = p, q = q)
# save.image('./data-raw/TGVAR.rdata')

# 输出阈值非线性检验和ADF检验
# xlsx::write.xlsx(Hansen99(countryrlt)$ThrTest,'E:/17_HuaDong/reserch/GVAR/RltPaper/ThresholdTest.xlsx',
#                    'Sheet5',append = T,row.names = FALSE)

# 计算负向脉冲冲击
impvar <-  'unctn'
impcnt <-  'USA'
responseCnt <- 'CHN'
resvar <- 'TFP'

irfdn <- GIRFTGVAR(countryrlt,regdata, endovar,TradeWeight,regime = 'dn', impvar = impvar,impcnt = impcnt, p = p, q = q)
irfup <- GIRFTGVAR(countryrlt,regdata, endovar,TradeWeight,regime = 'up', impvar = impvar,impcnt = impcnt, p = p, q = q)

# 画图
DrawIRF(dn = GetIRFPicdata(irfdn, cum = T, n = 24), up = GetIRFPicdata(irfup, cum = T, n = 24),
        responseCnt = responseCnt,resvar = resvar)

svfl <- paste('./data-raw/',impcnt,'_',EPUtype,impvar,'2',responseCnt,'_',indus,patnm,resvar,'_p',as.character(p),sep = '')
ggsave(paste(svfl,'.wmf',sep = ''))
save.image(paste(svfl,'.rdata',sep = ''))
# usethis::use_data(regdata, overwrite = TRUE)


